WKB approach and quantum corrections to classical dynamics in the Josephson 
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We apply a many-body WKB approach to determine the leading quantum corrections to the semi- 
classical dynamics of the Josephson model, describing interacting bosons able to tunnel between two 
localised states. The semiclassical dynamics is known to divide between regular oscillations, and 
self-trapped oscillations where the sign of the imbalance remains fixed. In both cases, the WKB 
wavefunctions are to be matched to Airy functions, yielding a modified Bohr-Sommerfeld quanti- 
sation condition. At the critical energy dividing normal and self-trapped oscillations, the WKB 
wavefunctions should instead be matched to parabolic cylinder functions, leading to a quantisation 
formula that is not just the Bohr-Sommerfeld formula, and recovering the known logarithmic quan- 
tum break times at this energy. This work thus provides another illustration of the usefulness of 
the WKB approach in certain many-body problems. 

PACS numbers: 03.75.Lm, 03.65.Sq, 03.75.Kk, 05.45.Mt, 
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I. INTRODUCTION 

The collective dynamics of a large number of interact- 
ing quantum systems can often be described semiclassi- 
cally, as mean-field approximations of the dynamics of 
such systems become more accurate with increasing sys- 
tem size. The great progress in trapping and manipu- 
lating cold atoms, and in producing strong coupling be- 
tween confined photon modes and matter degrees of free- 
dom have both led to an increasing variety of systems 
where it becomes possible to study many-body dynamics 
in isolated systems, and to investigate the extent to which 
semiclassical descriptions are applicable. In many cases, 
the classical dynamics can show or has shown interesting 
collective oscillations; examples include interconversion 
between fermionic atoms and molecules [lH3|, dynamical 
superradiance in coupled light-matter systems 0, @, op- 
tomechanical oscillations 6j , and the topic of this article, 
Josephson oscillations between atoms trapped in different 
wells 

Josephson oscillations [l2j of atoms between two 
trapped condensates are well described by the Hamil- 
tonian 



H = U{a)a - tfbf + J(a^b + tfa). 



(1) 



The semiclassical dynamics of this system have been 
studied in the context of the Lipkin-Meshkov-Glick 
Hamiltonian: Introducing the quantum operators S z = 
(tfb-rfa)/2,S x = (a^b + tfa)/2,S y = i(a?b-tfa)/2, the 
Hamiltonian becomes 



H = AUSt + 2JS X . 



(2) 



The semiclassical dynamics divide into two regimes: 
Small regular oscillations between the wells for low en- 
ergies and self-trapped oscillations, where the imbalance 
remains of a fixed sign for high energies (l3j . Specifically, 
if we replace the operators by c-numbcrs, the canonical 
equations of motion reduce to the single expression 



S z 



J Z N 2 



E 2 + {2UE~ J 2 ){2S Z ) 2 -U 2 {2S z f. (3) 



Here we have used conservation of energy and parti- 
cle number N 2 /4 = \S\ 2 . If J 2 N 2 - E 2 < and 
2UE — J 2 > then the region of allowed S z divides into 
two, excluding small population imbalances (small values 
of S z ). This means that for 2NU > J, then there are two 
classes of oscillations; those with energies E > N J are 
self-trapped oscillations 0, H, [U| where the sign of im- 
balance remains fixed, and for those with E < NJ the 
sign of the imbalance varies periodically. Both regimes 
have been observed experimentally 0, 03] • The Lipkin- 
Meshkov-Glick Hamiltonian, Eq. ^ , has been studied in 
a coherent state representation semiclassically [lH and 
with quantum corrections [15| . It has also been studied 
for small numbers of particles where quantum corrections 
may play a larger role [Tfj| . 

Instead of the spin representation, this article maps 
the Josephson problem onto a discrete Schrodinger equa- 
tion. As this is a one-dimensional problem, it can be 
approximately solved by the WKB approach 17]. Ex- 
panding for large system size N, the WKB wavefunc- 
tions can be found for the allowed and forbidden regions 
of Fock space. At the boundaries, these are matched to 
exact solutions of the Schrodinger equation, which often 
yields Airy functions. In such cases the spectrum is given 
by the usual Bohr-Sommerfeld quantisation condition on 
the action A(E n ) = J pdx = 2ir(n + 1/2) at a particu- 
lar energy level E n . Using this phase-space quantisation 
approach, the Josephson problem has been studied in 
pjl, [lH . One notable feature of the validity of semiclas- 
sics for the Josephson problem is that the quantum break 
times — the characteristic time at which quantum and 
classical dynamics start to differ — may grow only loga- 
rithmically, rather than algebraically, with when near 
the energy dividing self-trapped and non self-trapped os- 
cillations. This has been seen both by using a cumulant 
expansion [20l . l2ll ] and by phase-space quantisation [22j ■ 
The existence of an associated logarithmic divergence in 
the density of states has been discussed by Hooley and 
Quintanilla 23] for the related model of a single particle 
moving in an infinite tight-binding lattice with an overall 
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harmonic trap. 

The aim of this article is to apply the many-body WKB 
approach as discussed in Ref. 24] to the Josephson prob- 
lem of Eq. (JT|), and to use this to find the quantum break 
times for critical and non-critical levels. We solve the 
Schrodinger equation for large system size N to next-to- 
leading order. The WKB wavefunctions in the allowed 
and forbidden regions are then matched to appropriate 
special functions at the boundaries. For most energies, 
the appropriate special function is the Airy function, 
which then recovers the Bohr-Sommerfeld quantisation 
condition derived in [l8|, Ojj]. However, near the criti- 
cal energy separating self-trapped and non self-trapped 
oscillations, the relevant special functions are instead 
parabolic cylinder functions, and these lead to a differ- 
ent quantisation condition, that goes beyond the Bohr- 
Sommerfeld formula of [TH, [l9j]. Even away from the 
critical level, the full WKB analysis yields information 
beyond Bohr-Sommerfeld quantisation: The allowed re- 
gion divides into a regime where matching to a purely 
decaying Airy function is appropriate and one where the 
decaying solution is also highly oscillatory. While this 
behaviour modifies the nature of the wavefunctions, its 
effect on the quantisation condition is only a jump in the 
Maslov index. As such, it leaves the density of states 
unaffected, and so it was not crucial to previous work in 
the context of Bohr-Sommerfeld quantisation. 

The approach presented in this article has also been 
used in Ref. [24| to study the Tavis-Cummings model [25[ 
describing e.g. the dynamics of two-level systems cou- 
pled to photons. The many-body WKB approach gives 
quantum dynamics for this model very similar to the clas- 
sical dynamics for most energy levels and initial condi- 
tions. However, there are critical levels, for which the 
matching of the WKB wavefunctions becomes more com- 
plex and gives rise to anharmonicity scaling as the loga- 
rithm of system size [2(|. In both the Tavis-Cummings 
model, and in the Josephson model, these critical levels 
show logarithmic quantum break times, analogous to the 
ln(7i) quantum break times found near unstable classical 
states in single particle quantum problems [13] . A fur- 
ther class of problems in which this many-body WKB 
approach may be of use concerns variations of problems 
such as the Josephson problem, and the Tavis-Cummings 
model, in which parameters are varied as a function of 
time to give many-body generalisations of Landau-Zener 
problems[28-33]. Many treatments of this problem have 
been effectively semiclassical, and the WKB approach 
may provide a method to determine whether large quan- 
tum corrections can ever arise as a result of transitions 
to and from critical levels. An approach along these lines 
has been explored by Witthaut et al. 
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The remainder of this article is arranged as follows. 
Section |H] maps the many-body Hamiltonian to a form 
amenable to solution by the WKB method, and then pro- 
vides the ingredients necessary to determine this solution. 
These ingredients are the WKB wavefunctions in the clas- 
sically allowed and forbidden regions, given in Sec. Ill Al 



and the connection formula holding at the boundaries be- 
tween these regions, given in Sec. Ill B[ Section IIIII then 
combines these ingredients to give the resulting quanti- 
sation conditions, which are specified in three separate 
ranges of energies, according to whether one is above, 
near, or below the critical energy level (the energy level 
dividing self-trapped and non-self-trapped oscillations). 
From these quantisation conditions, section IIVI then ex- 
tracts the scaling of the quantum break time as a function 
of system size, and summarises the results found. 



II. DISCRETE WKB APPROACH TO 
JOSEPHSON EQUATION 

To use the WKB approach, we need to produce a dis- 
crete Schrodinger equation, by writing the state of the 
system in terms of the total number of particles N, and 
the imbalance n: 



N 



-N 



a UN+n)/2 b UN-n)/2 

^[{N + n)/2]![(JV -n)/2}\ 



\Q) (4) 



where n is restricted to the same odd/even parity as N. 
Acting on this state with the many-body Hamiltonian, 
and looking for eigenstates with energy E, yields a dis- 
crete Schrodinger equation: 



(E - Un 2 )ij; n 



J r 

2 



y/(N + n)(N -n + 2)i, n _ 2 
+y/(N + n + 2)(N- n)i/j n+2 ] . (5) 



Writing E = eJN, U — Ju/N, z = n/N, one can sepa- 
rate the system size dependence from the other parame- 
ter dependence and write: 



N 



N 




1-z) 1 + z 



(6) 



In these units, self-trapped states exist only if u > 1/2 
and e > 1. 



A. WKB wavefunction 

For a given energy, one may divide the range of — 1 < 
z < 1 into classically allowed and forbidden regions. 
These are distinguished by oscillating versus decaying 
wavefunctions, and correspond directly to the regions of 
S z = Nz in Eq. ([3]) which the classical dynamics explores. 
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1. Allowed region 
In the allowed region, the WKB ansatz has the form: 



ij(z) = b(z) 



Q e KNW +W x ) + c _ e -i(NW +W 1 ) 



(7) 



where Wq and W\ are the z-dependent phase terms at 
leading and next-to-leading order. By substituting this 
into Eq. Q and identifying real and imaginary terms of 
the same order in 1/N, one finds the definitions: 



same parity as N exist, this corresponds to terms having 
alternating signs. After this transformation, Eq. (|51- [TU)) 
are modified by replacing (e — uz 2 ) — > — (e — uz 2 ). The re- 
gions where this transformation is necessary are indicated 
in Fig. [TJby hatching. We will use the notation Wq, W\ 
for the phases calculated with this additional minus sign. 
When part of the allowed region does require this trans- 
formation, and part does not (i.e. for < e < u), it is 
necessary to consider carefully the connection between 
these regions; this is discussed further in Sec. Ill B 31 



cos(2Wq(z)) = 



VT^ 2 
b(z) = [l-z 2 - (e-uz 2 ) 2 ] 

9 

W[(z) 



-1/4 



2(1- z 2 )y/l- z 2 - (e 



uz* 



(8) 
(9) 
(10) 



(Note that the sign written for W[ assumes that the 
solution of Eq. © is taken such that sin(2WQ) > 0). 
For this solution to be valid, it is clearly necessary that 
\e — uz 2 \ < \Jl — z 2 , which defines the classically allowed 
region shown in Fig. [1] 
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u=1.2 



FIG. 1: (Color online) Boundaries between allowed and for- 
bidden regions, and between regions which needs extra fac- 
tor of (-l) n/2 . Solid shaded region is classically forbidden, 
hatched region (e < uz 2 ) requires extra factor of ( — l) n / 2 to 
allow matching to connection formula at boundary. Plotted 
for u = 1.2 

Within the classically allowed region, there is also a 
division between e > uz 2 for which < 2Wq < ir/2, and 
e < uz 2 for which tt/2 < 2Wq < tt. While this distinction 
is unimportant within the allowed region, it will later be 
necessary to distinguish these regions, so that the connec- 
tion formulae at the boundary of the allowed region can 
be written in terms of smooth wavef unctions. If e > uz 2 , 
then at the boundary Wq — > 0, so the wavefunction is 
smooth. For e < uz 2 , Wq — > tt/2, which indicates that 
NWq = iNxw/2, which is rapidly varying. This rapid 
variation can be removed in such cases by instead defin- 
ing i[) n — ¥ ipn(— I)™/ 2 - Since only values of n with the 



2. Forbidden region 

In the forbidden region, the wavefunction exponen- 
tially decays, so the WKB ansatz becomes: 

i/)(z) = b(z) \c +e ( Nn ° +n ^ + C_e^ Nil °+M . (11) 

In this case, there is no strict distinction between the 
terms to be incorporated in b(z) and fii(z) — both de- 
scribe the real part at order 1/N. However, it is con- 
venient for the connection formulae to keep b(z) as in 
Eq. ([9|), after which J7o,i can be identified by powers of 
1/N as: 



cosh(2fi[,) = 

n[ = 



— (e — uz 2 ) 



2(1 - z 2 )y/(e-uz 2 ) 2 - 1 + z 2 



(12) 
(13) 



In this case, the sign of assumes Sl > 0. In the 
forbidden region at small z one always has e > uz 2 so no 
alternating sign factors are needed, but in the forbidden 
region at large z, a factor (— l)™/ 2 may be needed so 
that cosh(2f2g) > 0. The distinction between exponential 
decay, and exponential decay with alternating signs is 
clearly visible in the exact wavefunctions shown in Fig. [21 
In the forbidden region at small z, one will have |C_| = 
|C+| (assuming that f2 is measured from z = 0), since the 
wavefunctions will be either odd or even functions of z. 
In the region at large z, only one of C'± will be non-zero, 
describing exponential decay for z — > =poo. 



3. Critical level 

At the critical level e — 1, there is a bifurcation, and 
the WKB wavefunction in the allowed region must be 
matched to a different special function for z ~ 0. For 
this purpose it is convenient to write e = 1 + X/N, where 
A is small compared to TV, and to rewrite the WKB form 
for the allowed regime making use of this — because the 
deviation of the energy from 1 is now considered as being 
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level, the eigenvalue spacing will depend logarithmically 
on the system size, which in turn leads to logarithmic 
dependence of the quantum break time on system size. 



B. Connection formulae 

To connect the WKB ansatz valid in the classically 
allowed and classically forbidden regions, one needs the 
solution valid near the boundary. If the wavefunction 
is smooth near this point, one may expand Eq. ([5]) di- 
rectly in powers of 1 /AT to give a differential equation for 
the wavefunction. (Note that in the WKB ansatz, there 
was no such assumption of smoothness, since we did not 
require Wq to be small.) To ensure the wavefunction 
is smooth, a factor (— l)™/ 2 may be needed according to 
whether the boundary is in the hatched on unhatched re- 
gion of Fig. [3 as discussed above. Both these cases can be 
considered together, by replacing (e — uz 2 ) — » |e — uz 2 \ to 
take account of the two possible signs. With this change, 
expanding Eq. (|BJ) to order 1/N 2 one finds: 



uz 2 \ip 



Vl-z 2 + 



N 2 



zip' 



O (N- 



(17) 



When considering the leading and next to leading or- 
der behaviour (i.e. classical behaviour plus leading order 
quantum corrections), the term of order jV -2 term in- 
volving ip may be neglected, as there are larger terms in- 
volving ip, whereas the terms involving ip' and ip" should 
be kept. Rewriting the above equation one has: 



= i\>"- 



N 



[2(1 



iV 2 
2 



Vi^ 2 



(18) 



FIG. 2: (Color online) Wavefunctions at various characteristic 
energies, showing behaviour in allowed and forbidden regions. 
Background shading is as in Figff] Plotted for u — 1.2 



of order 1/N, its effect is relegated from Wq to W\, 
1 — uz 2 



cos(2Wq(z)) = 



VT^z 2 

b(z)=[z 2 (2u-l-u 2 z 2 )Y 1,A 
1 - A + (A — u)z 2 



W[{z) = 



2(1- z 2 )\zW2u~l-u 2 z 2 



(14) 

(15) 
(16) 



One can see from Eq. (116)) that Wi will be logarithmically 
divergent as z — > 0. The form of ip(z) for z ~ will give 
the natural cutoff to this logarithm, which will be found 
to depend on the system size N, so exactly at this critical 



1. Regular boundaries — Airy functions 

Away from the critical level e = 1, the boundary 
between the allowed and forbidden region is given by 
|e — uz 2 \ = s/l — z 2 , which has the outer (inner) solu- 
tions: 



1 



z 2 . = 

°' 1 2u 2 



2ue - 1 ± y/l - Aue + Au 2 



(19) 



Near these boundaries, we want to write Eq. (|18[) in terms 
of distance C, from the boundary, z — ±z 0) i + C- Since 
Eq. (|18j) is even under z —> —z, it is clear that the equa- 
tion at z — z ,i + C an d z — ~ z o.i + C are related by 
( — > — £. Therefore, considering the former of these, we 
may write Eq. (fl8|) as: 



= V"-a o , i V / + [/3 o ,i + /o,i(C)]V' 



(20) 
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with a o4 = z 0>i /(l - (),&,< = JV/[2(1 - and: 



/o,i(0 



TV 2 



\e - uzli - 2uz ,i(\ 



2z G v£ 



2u{e - uz 2 Q ,) - 1 



.7V 2 2o ,CVl-4ue 



del' 



7o,iC (21) 



In the above, we have assumed that £ is small so the sign 
of e — mz 2 does not change, and the last line has used the 
form of z 2 i in Eq. (|19p , with the =F signs corresponding to 
the outer (inner) boundary. With /(C) = 7Cj the solution 
Eq. (|20[) can be written using Airy functions: 



-«C/2 



C o A(-7 1/3 + C bJ Bz(- 7 1/3 e)l (22) 



with £ = £ + /3/7 — a 2 /A^j. For the outer boundaries we 
match to a decaying solution, so Cb = 0. Since 7 Q is neg- 
ative, the solutions is oscillatory for £ < and decaying 
for £ > as expected. For the inner boundaries, both 
exponentially decaying and growing parts are required. 
Since 7, is positive, the solutions are oscillatory for £ > 
and grow/decay for £ < 0. 



2. Critical boundary — Parabolic Cylinder functions 

Near e = 1, the inner boundary becomes an extremum 
at z = 0, and the form of Eq. (|18l) is different to that in 
Sec. lII BTl To study energies near this level, we will write 
e = 1 + X/N as in Sec. Ill A 31 Since e = 1, z = never 
requires a factor (— l)™/ 2 in the wavefunction, Eq. (TI5)) 
can always be written near this point as: 



= i/}"- zip' + 



N N 2 
_(1-A) + — (2u-l)*' 



^ (23) 



To solve this equation, it is convenient to define: 



fi = y/2u — 1, 



X = 



1 - A 
2/x 



(24) 



After removing a Gaussian factor ip = e z2//4 /, this can 
be recognised as Weber's equation[3J|; in terms of £ = 
'iV/x one has = /" + («x— £ 2 /4)/ with solutions 



re 



-i7r/4 



in terms of parabolic cylinder functions. The general 
solution can be written as: 



- ^ 2 /4 



oAx-i/2 (e~ W4 V^) 



(25) 



The Gaussian prefactor only contributes at order 1/N, 
so in matching the asymptotics of NWq, N&Iq to ijj, this 



prefactor can be dropped. It is clear that for z — > —z, this 
expression changes as a -H> /?. The asymptotic expansion 
of this expression for large x given in Ref. [34[ can be 
written: 



aexp(^ — +*- 



+ /3exp 



x exp 



■P 



N/iz 2 



In 



3ttx _ 
4 

Zy/WJl 



.3tt 



2tt 



exp 



r[i/2 - i X ] 

Nfiz 2 



x exp — 1 



7TX 7T 

' 4 *8 
1 

-ix- 5 



In 



(26) 



As noted in Sec. Ill A 31 the phase has a logarithmic di- 
vergence at small z, and the form of the wavefunction 
given here provides the cutoff for the logarithm, 1 / y/Nfi, 
which depends on system size. 



3. Connection within allowed region 

In addition to connection formulae at the boundaries of 
allowed and forbidden, for energies in the range 1 < e < 
u, it is necessary to connect solutions with and without 
the extra factor of (— l)"/ 2 in the middle of the allowed 
region. 

In such cases, we may write: 



(_l)n/a 



C+e iW 



C_e- lW 



C + e iW 



C-e- iW 



\z\ < z s 
\z\ > z s ' 



(27) 



where z s = \feju is the point at which the sign of e — uz 2 
changes. The question is how to relate C± to C± at 
the boundary. To fully define the the phases W, W one 
must specify the limits of integration. It is convenient to 
choose: 



W{z)= I dz[NWo + W{], W{z)= I dz[NWo + W{] 

(28) 

and to then define the phases at z s as AWi n = 
W(z s ), AW out = —W(z s ). To determine the connection, 
we focus on the phase due to Wq] expanding for z = z s +£ 
one has: 



cos(2W^) = 



-2uyi£ 
y/u — e 



cos(2W/o), (29) 



and integrating one has: 
W(z s + £) = +AW in + N 

W(z s + £) = -AW ont + N 



77 » 
4 C + 



c 2 



VT~^l 2 



(30) 
(31) 
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To match the £ dependence near £ = one can write the 
prefactor (-1)™/ 2 = exp[inN(z s + ()/2] which gives the 
matching conditions: 



C_ = C + exp 



JAW;, 



■AT?r 



(32) 



and the complex conjugate for the other pair. 



III. BOUNDARY MATCHING AND 
QUANTISATION CONDITIONS 

By connecting the wavefunctions in the allowed and 
forbidden regions to the Airy or parabolic cylinder func- 
tions at the boundaries, one can derive conditions relat- 
ing the coefficients C± in these various regions. In order 
to simultaneously satisfy all the equations, it is neces- 
sary to have certain conditions on the WKB phases; this 
provides quantisation conditions on the allowed energies. 
This section will derive these quantisation conditions in 
the three cases e < 1, e > 1 and e = 1 + X/N. In all 
cases the quantisation condition can be reduced to a form 
/(e) = mir, where m is an integer. Figure[3]illustrates the 
accuracy of these WKB quantisation methods by plot- 
ting min m (/(e„) — mir), where e n is an exact eigenvalue 
— i.e. determining the extent to which the exact eigen- 
values obey the WKB quantisation condition. 
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FIG. 3: (Color online) Accuracy of WKB-derived quantisa- 
tion condition. Each line corresponds to difference of appro- 
priately defined phase difference, and nearest integer multiple 
of 7T, using quantisation conditions appropriate to energies 
below, near and at the critical level. Plotted for u — 1.2 



A. Below critical level 

We consider separately the cases e < and < e < 1. 
For e < 0, the quantisation condition is simple as the 
factor (— l)™/ 2 is needed throughout the allowed region. 



Let us define the limits of integration for W(z) such that: 

W(z) = [ dz[NWo + W{], AW = W(z Q ). (33) 
Jo 

Then, by expanding the equation for Wq (i.e. Eq. (|T4|) 
with an extra minus sign) near z = zo + Ci an d comparing 
the result to the last line of Eq. (|2"Tj) one can show that: 



/\2 



2u(e — uz"l) — 1 



1-zi 



ZoC = 1- 



27o 
' N 2 



C (34) 



Remembering that 7„ < 0, integrating the expression for 
Wq gives: 

inz +c) = AVF-^i 7o r /2 ici 3/2 , (c<o). (35) 

In order to match ip = (— 1)™/ 2 [C' + exp(iW) + 
C- exp(— iW)] to the expansion of the Airy function [34| 
appearing in Eq. 



^cx^l^l 1 ^-!-...)] 



l7o 



l/2|£|3/2 



(36) 



gives the condition (C+ / C-) exp(2iAW) — i. Repeating 
this procedure at z — —z a + one has instead 



^z + C)=-A^+l| 7o | 1 / 2 | C |3/ 2; 



(C > 0), (37) 



and similar matching gives (C+/C_) exp(— 2iAW) = —i. 
Combining these, the quantisation condition is: 



2AW 



TV 

2' 



(38) 



which is just Bohr-Sommerfeld quantisation as expected 
for this simplest case. 

In the range < e < 1, there are both regions with 
and without factors (— l)™' 2 , however since both bound- 
aries have this additional factor, it is possible to get away 
with Eq. (f38|). using the integral in Eq. (|33|) even when 
cos(2VFq) < 0. To show that this is valid, we will dis- 
cuss the result of taking the connection at z = z s into 
account explicitly. Using the results of Sec. Ill B 31 means 
that in the connection formulae at the outer boundary, 
one should replace 



Q e iAW _^ (j e i[AlV out -AW in +7rjVz s /2] 



(39) 



and for C_ the condition is similar, but complex con- 
jugate. Here, the expression W m in Eq. (fBT)]) should be 
understood as taking zi — > in the expressions following 
Eq. (f2"5|) , since there is no inner boundary for e < 1. 
With this replacement, the quantisation condition be- 
comes: 



TO7T = 2 



AW n - 



AW h - 



irN 



7T 

2' 



(40) 
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To relate the integrals AW out ,AW in to the AW of 
Eq. (j3"3")l we should note two features. Firstly there is 
a term N(ir/2)z s in the difference of definitions due to 
the inverse cosine, which is compensated by the last term 
in brackets in Eq. (|40[) . Other than this, one may see 
that the various sign factors in front of AW[ n , in the 
relative definition of AWi n vs AWout, and in the or- 
der of limits of integration in Eq. ((25)) are such that 
the remaining parts of the integrals all match so that 
AW in - nNz s /2 - AW ont = AW. 



B. Above critical level 



boundary, one must match both Airy functions, which 
gives the matching condition: 



(ggj ~ iC R , a )e+"/ 4 = R+ 
(C R , b + iC R , a )e-™/* R_ 



(44) 



where C R ^ a ,C R ^ are the coefficients of Airy functions, 
as in Eq. (|22|) at the right hand (+Zi) boundary. These 
should also be matched to the coefficients of the grow- 
ing/decaying terms in the forbidden region. Defining 
Aflp = 2Clp[zi) — Qp(zi) — zi), one may expand 

n F ( Zl + () = (Afif/2) - (2/3)| 7l | 1 / 2 |C| 3/2 for C < 0. 
Matching this to Airy functions yields: 



Above the critical level, there can either be three or 
five separate regions to consider: a forbidden region 
in the centre and either single allowed regions to the 
left and right if e > it, or pairs of allowed regions if 
e < u. We will label these regions, from left to right as 
L, L, F, R, R, with the L,R regions having z s < \z\ < z a 
when e < u, and vanishing otherwise. These labels will be 
used both for the coefficients L±,R±, and for the phases 
Wi,Wl,Qf,W r ,W r . Focusing on the right hand side 
(z > 0), it is convenient to define W R , W R with the same 
limits as W,W in Eq. (|28[) . and to define ilp(z) as the 
integral from to z. 

Consider first the case e < u. With the above def- 
initions the outer boundary condition simply becomes 
= i, as the phase in this region is defined so 
it vanishes at z Q . Using Eq. (|32|) . this translates to a 
condition: 



R+ 

R- 



exp i2 



AWin ~ AW„ 



nN 



i, (41) 



where the tt/2 term on the left hand side is associated 
with an overall minus sign from i —> 1/i on the right 
hand side. 

In the case e > u we need only W(z), which we have al- 
ready defined as the integral from Zj to z. In this case the 
phase matching condition at the outer boundary implies: 

R f z ° 

— - exp (i2AW) = i, AW = / dz[NW^ + W[}. (42) 

Both Eq. (T4*l"j) and Eq. (|4"2")) can be combined by the state- 
ment exp(2iAW eff ) = i, with 



AW t 



cff 



AW in - AW Q 
AW 



EK Z -L 2L e < ii 

e > u 



The equivalent analysis at the leftmost boundary gives 
(L + /L_)exp(-2iAWcff) = -i. 

The boundary at Zi always involves R± rather than 
R±, and our choice of limits means that the phase is 
always measured from this boundary. Following the same 
logic as led to Eq. (|57|) . and noting that % is positive, 
one has W(z t + C) = (2/3)|7,| 1 / 2 |C| 3/2 for C > 0. At this 



C 



F _ e -AQ F /2 



C R , a F+e+ A0 ^/ 2 



R,b 



(45) 



Equations (|44|) and (|45|) can be combined to eliminate 
the coefficients of Airy functions giving: 



F _ e -AQ F /2 _ lF+e +A(l F /2 ^ R 
F _ e -AQ F /2 _|_ ip +e +ACl F /2 ~ 



+ -itt/2 _ e -2iAW off 



Analysis at the left hand boundary is very similar, 
p +e -An F /2 _ iF _ e+ An F /2 ^ L _ _^ 



(46) 



= , -/2 = e -2iAW«H 



p +e -AQ, F /2 + i p_ e +An F /2 p 

_ (47) 
Eliminating F± from these equations then yields the 
quantisation condition: 



1 + cos(2AW cff ) = _2AO f 
1 - cos(2AW eff ) 6 



(48) 



Since the quantity Aflp grows linearly with N, the quan- 
tum tunnelling between allowed regions is hugely sup- 
pressed. Taking this tunnelling into account at leading 
order yields cos(2AW c ff) ~ — 1 + 2 exp(— 2Ailp), so that: 



(49) 



Differentiating this gives an expression for the tun- 
nelling induced energy splitting (the inverse of quan- 
tum tunnelling time for self trapped states): <5e sp ntting = 
2e~ AnF \dAW cff /de\- 1 . 



C. Near critical level 

Near the critical level, we replace the matching to 
Airy functions and the forbidden region by matching to 
Eq. (|26[) . As in the previous section, we either have four 
allowed regions or two depending on whether e = 1 is 
less than or greater than u. Just as in that case, we may 
write the boundary conditions at the outer boundaries 
as: 



R 

R- 



+ e +2iAW Bt{ 



L 



± e -2iAW Bli 



(50) 
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with A W e s given in Eq. (|49|) , and Zi — > so the phase fac- 
tors vanishing at the inner boundary. However, because 
the expression for W± given in Eq. (Tpo)) diverges loga- 
rithmically at the critical level, it is necessary to slightly 
modify the definition of AW c g, so the lower limit of the 
integral for W[ is taken to be z cu toff rather than zero, to 
avoid the logarithmic divergence. 

Expanding Eq. (fLT| . (fPo)) near z = gives: 



l-2(^) 2 = l- l -{2u-l)z\ W{ = 



1-X 



2\z\y/2u- 1 



(51) 



Using the definitions of x, M given in Eq. (|24|) . these ex- 
pressions can be integrated to give: 



W R{L) (z) = ±^± X \n 



•^cutoff 



(52) 



for z > (z < 0) respectively. If we choose z cuto s = 
1/y/Npi, the coefficients R±, L± can then be matched to 
Eq. for z > (z < 0). This gives: 



/X 71 " 7r 
i?+ = aexp(^ — + i~ 



3xtt .3tt 



R- = 



r[| - ix] 

,e X p(f 



exp 



/ X 71 " . ■ ^ \ 



a exp 



3xtt 
4 



.3tt 



ay 27r 



exp i 



(53) 
(54) 
(55) 
(56) 



By making use of the Weierstrass identity [34| to write 
r[| — «x]r[^ + ix] = 7r/ cosh(x7r), these expressions can 
be combined to give the quantisation condition: 



Arg 



'2it 



r[i/2 + *x] 



e 2 ' 4 ^ + e 



IV. CONCLUSIONS AND SCALING OF 
QUANTUM BREAK TIME 

We have shown how semiclassical quantisation formu- 
lae emerge from a many-body WKB approach for the 
Josephson problem. We now use this to extract the 
quantum break time; this time is given by the inverse 
of the anharmonicity of the energy level spacing. At low 
and high energies Bohr-Sommcrfeld quantisation applies 
with small 1/N quantum corrections. In between these 
two regimes, a critical level with large quantum correc- 
tions is seen, where larger deviations from semiclassics 
occur, leading to quantum break times that scale with 



the logarithm of the system size. From the quantisation 
formulae in the various regimes, we can straightforwardly 
extract the anharmonicity of the spectrum. For perfectly 
regular spacing of energy levels, the dynamics would be 
periodic, matching the semiclassical dynamics. For the 
critical level, the energy appears in the quantisation con- 
dition only through A = N(e—1) as a prefactor appearing 
in Wi, and through x = (l~ A)/2/i in the gamma func- 
tion in Eq. (j57|) . The former of these contributions can 
be written as: 



^"AWeff 

aX 



=fi/v 



dz 



_2zy r , 



2/x 



Zcutoff =\/y/Wji 

(58) 

This logarithmic dependence means that the level 
spacing near this critical level is Se = 5X/N ~ 
Tr/[NdAW cS /dX] - l/[AHn(A0], by making use of 
Eq. (|58]l . If this logarithm were large, then the solutions 
of Eq. (|57|) would require A to be regularly spaced. How- 
ever, because of the A dependence in F[l/2 + ix], there 
will be some anharmonicity. Since the only parameter 
controlling the anharmonicity is dAW c g/dX, the anhar- 
monicity of level spacing SSe = e n +i + — 2e„ is given 
by SSe = 55X/N ~ l/[iVln(A^) 2 ]. In contrast, away from 
the critical level, level spacing is 1 /N and anharmonicity 
1/N 2 . Thus, the quantum break time scales logarithmi- 
cally with system size near the critical level and linearly 
elsewhere. This is exactly the conclusion found by con- 
sidering next to leading order corrections to the quan tum 
dynamics via a cumulant expansion in Ref. [U Ml and 
consistent with Ref. [22l |. 

In conclusion, we have shown how the WKB approach 
for many-body systems, discussed in 0, Hfl can be ap- 
plied to the dynamics of the Josephson problem, illus- 
trating that it can be applied to another paradigmatic 
problem of collective quantum dynamics. This many- 
body WKB approach is particularly useful in cases where 
critical energy levels exist. In such cases, the semiclas- 
sical description may be inadequate even for mesoscopic 
systems with up to ~ 10 6 particles, yet such numbers of 
particles make numerical approaches to the full quantum 
dynamics very expensive. As such, it provides an ideal 
tool to identify cases where quantum dynamical effects 
survive in mesocopic systems. 
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